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Section  1  Introduction 

1.1  Research  Summary 

This  document  reports  on  the  work  carried  out  by  Honeywell  Laboratories  from  1999  to  2001  under 
the  AFOSR  contract  F49620-99-C-0009  entitled  “Collective  Management  of  Satellite  Clusters.” 

During  this  contract  we  carried  out  work  in  three  areas:  orbital  analysis,  cluster  reconfiguration, 
and  cluster  control.  These  areas  correspond  to  the  top  two  layers  of  the  proposed  hierarchy  for 
cluster  management  shown  in  Figure  1.1.  We  did  not  address  problems  relevant  to  the  bottom 
layer,  which  is  concerned  with  attitude  control  and  regulation.  It  is  important  to  note,  however, 
that  this  hierarchical  decomposition  is  to  some  degree  arbitrary  and  that  it  leads  to  a  certain  loss  of 
performance  in  exchange  for  design  simplicity  and  robustness. 

Method  Signal  Layer  Signal 

Human  in  the  Loop  Mission  Command 


Physical  Layer 


Figure  1.1:  Control  hierarchy  for  cluster  management 

We  approached  the  subject  of  orbital  analysis  for  cluster  maintenance  from  two  points  of  view. 
In  the  work  reported  in  Section  2,  we  define  an  approximate  separation  function,  relatively  simple 
to  compute,  to  help  predict  the  short-term  behavior  of  the  constellation  and  to  aid  in  tasks  such  as 
collision  avoidance.  We  do  further  orbital  analysis  in  Section  6,  where  we  show  how  direct  opti¬ 
mization  techniques  can  be  used  to  generate  orbital  paths  that  meet  mission  objective  requirements 
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1.2.  Personnel 


(cluster  shape)  and  require  minimal  fuel  usage  for  their  maintenance. 

In  Section  3  we  show  how  market-oriented  programming  (MOP)  methods  can  be  applied  to 
the  problem  of  cluster  reconfiguration.  Cluster  reconfiguration  is  necessary  because  maintaining 
operational  formation  for  long  periods  of  time  can  be  prohitively  expensive,  It  is  UiUS  necessary 
to  move  the  constellation  from  parking  orbit  into  mission  orbit  and  back  with  minimal  u^e  of  fuel. 
Although  suboptimal,  MOP  is  extremely  simple  to  implement  and  yerify,  even  when  the  cluster 
has  a  changing  number  of  components.  In  this  section  we  also  analyze  the  robustness  of  distributed 
management  algorithms  to  different  types  of  failure  modes  to  guarantee  their  safety. 

Finally,  we  cover  the  issues  relevant  to  cluster  guidance  in  Sections  4  and  5.  In  these  sections 
we  present  an  approach  to  cluster  guidance  that  achieves  very  low  levels  of  fuel  usage  while  being 
easily  ammenable  to  all  the  operational  restrictions  of  the  Techsat21  mission.  In  section  5  we 
discuss  a  detailed  study  of  the  effects  of  different  factors,  such  as  sensor  noise  or  constellation 
shape,  on  the  total  fuel  requirements  of  the  cluster.  For  this  study,  we  use  the  parameters  of  the 
Techsat21  demonstrator  mission,  planned  for  2003,  to  analyze  its  feasibility  and  further  direct  its 
development. 


1.2  Personnel 

The  following  scientists  have  been  partially  supported  by  this  contract:  Dr.  Jorge  E.  Tiemo,  Dr. 
Blaise  Morton,  Dr.  John  Samson,  Nicholas  Weininger,  Todd  Carpenter,  Scott  Snyder,  and  Lynn 
Gravatt. 


1.3  Publications 

“Collective  Management  of  Satellite  Clusters”,  in  Proceedings  of  the  1999  AIAA  Conference  in 
Guidance  and  Control. 

“Control  of  LEO  Satellite  Clusters”,  in  Proceedings  of  the  2000  IEEE  Conference  on  Decision  and 
Control. 
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Section  2  Orbital  Analysis  and  Satellite  Cluster 
Maintenance 


Our  approach  to  orbital  analysis  is  based  on  instantaneous  Kepler  orbits.  The  true  paths  followed 
by  satellites  near  the  Earth  are  not  Kepler  ellipses,  owing  to  forces  associated  with  the  non-point- 
mass  nature  of  the  geopotential  (e.g.,  the  J2  term),  the  moon’s  gravity,  atmospheric  drag,  and  other 
factors.  The  point-mass  geopotential  dominates  these  other  effects,  usually  by  several  orders  of 
magnitude,  so  for  short  periods  of  time  the  Kepler-ellipse  approximation  is  a  reasonable  choice 
for  the  nominal  motion,  about  which  perturbations  may  be  superimposed.  In  fact,  the  Kepler- 
ellipse  approximation  may  become  extremely  accurate  when  used  to  forecast  the  relative  motion 
of  identical  spacecraft  flying  in  close  proximity  (for  times  on  the  order  of  one  orbital  period), 
because  the  perturbative  forces  tend  to  affect  all  the  spacecraft  in  the  same  way.  On  the  time  scale 
of  a  single  orbital  period,  the  perturbations  will  cause  the  Kepler  orbits  of  the  various  satellites  to 
drift  together. 

When  the  time  scale  is  enlarged  to,  say,l(X)  orbital  periods,  however,  we  expect  to  observe 
divergences  in  the  parameters  of  the  Kepler  ellipses  of  the  various  satellites.  When  the  divergences 
become  sufficiently  large,  it  will  be  necessary  to  apply  corrective  thrust  to  one  or  more  spacecraft 
to  bring  the  configuration  back  together.  For  example,  if  the  energies  of  the  orbits  drift  apart,  there 
will  be  no  way  to  keep  the  flock  together  because  the  orbital  period  is  a  function  of  the  energy. 
Thus,  a  part  of  the  control  strategy  must  be  to  keep  the  energies  of  the  orbits  close  together. 

To  formulate  the  complete  control  strategy  the  following  issues  must  be  addressed; 

1 .  How  do  we  monitor  the  state  of  the  satellite  system? 

2.  What  are  the  desirable  states  for  the  system  to  be  in? 

3.  How  do  we  determine  when  a  thrust  maneuver  is  necessary? 

4.  When  necessary,  how  do  we  select  the  appropriate  maneuver? 

One  step  toward  answering  these  questions  is  to  develop  and  analyze  a  separation  function  for 
pairs  of  satellites.  The  separation  function  tells  how  close  together  and  how  far  apart  each  pair  of 
satellites  is  expected  to  get  during  the  next  period  (using  the  Kepler-ellipse  approximation).  The 
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2.1.  Notation  and  Coordinates 


values  of  the  separation  function  for  all  pairs  of  satellites  will  provide  one  way  to  monitor  the 
configuration.  It  may  be  necessary  to  consider  other  functions  as  well,  but  certainly  we  need  to 
monitor  bounds  on  intersatellite  distances  to  keep  the  system  healthy. 

To  decide  when  maneuvers  are  necessary,  we  will  need  to  determine  acceptable  ranges  for  the 
separation  functions.  For  example,  it  might  be  the  case  that  we  want  to  keep  the  distance  between 
satellites  at  least  10  m  but  no  more  than  1000  m.  One  question  to  address  in  passing  is  whether  a 
specified  range  is  even  achievable.  In  any  case,  the  basic  idea  is  to  mandate  a  maneuver  whenever 
one  pair  of  satellites  is  outside  the  acceptable  range  for  the  separation  function. 


2.1  Notation  and  Coordinates 

At  each  instant  of  time,  a  satellite’s  position  and  velocity  vectors  are  associated  with  the  parameters 
of  a  Kepler  orbit  about  the  Earth’s  center  of  mass  consistent  with  Newton’s  law.  We  will  be 
working  with  orbiting  satellites  for  which  the  position  vector  x  and  the  velocity  vector  v  =  ^ 
are  linearly  independent,  so  the  angular  momentum  h  (per  kilogram)  defined  by  h  =  x  x  v  is 
a  nonzero  vector.  If  the  point-mass  geopotential  were  the  only  force  acting  on  the  satellites,  then 
each  satellite’s  h  would  be  time  invariant.  The  instantaneous  Kepler  motion  of  each  satellite  is 
confined  to  the  plane  through  the  linear  subspace  normal  to  its  h  vector. 

2.1.1  Velocity  Circle 

Pick  a  coordinate  system  x  =  {x,  y,  z)  such  that  h  points  along  the  positive  z-axis.  Then  from 
x(t)  X  v(t)  =  h  =  (0, 0,  h)  we  see  that 

x(f)  =  (r(t)cos(O(t)),r(t)sm(0(t)),O) 

Differentiating  this  last  equation  with  respect  to  t  and  taking  the  cross  product  with  x(t),  we  find 

(0,0,/t)  =  h  =  X  X  ^  =  (0,0, 

So  in  particular, 

f 

Now,  from  Newton’s  equation  ^  where  k  =  GMg  and  setting  i?  =  A://i,  we  obtain 

after  integration: 

V  =  i?(-sin(0),cos(0),O) -fc 


Use  or  disclosure  of  data  contained  on  this  sheet  is  subject  to 
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2.2.  Characterization  of  Orbits  with  Fixed  Eneigy  and  Eccentricity 


where  c  =  (ci,  C2, 0)  is  an  integration  constant.  Interpreting  this  last  equation,  we  see  that  the 
velocity  vector  v  moves  along  a  circle  centered  at  c,  with  radius  R  =  k/h,  and  lies  in  the  plane 
through  the  origin  which  is  orthogonal  to  h.  Also,  note  that  v  —  c  is  always  orthogonal  to  the 
position  vector  X  =  r(cos(0),5in(^),  0).  (See  [3]  for  more  details.)  Specifically, 


2.1.2  Eccentricity,  True  Anomaly,  and  Energy 

Defining  c  =  l|cl|/i?,  choose  coordinates  (uju)  for  the  plane  of  the  velocity  circle  so  that  the 
center  of  the  circle  lies  on  the  positive  u-axis.  Then  c  is  called  the  eccentricity  and 

V  =  i?(— sm(0),€  +  cos(^),0) 

Substituting  into  h  =  x  x  v,  we  obtain 

h  =  ri?(l  +  ecos{0)) 

and 

r  =  A/(l  +  €cos{0))where  A  =  h^/k 

The  parameter  6  defined  as  above  is  commonly  called  the  true  anomaly  of  the  satellite.  Perigee 
(r  =  rjnin)  occurs  when  ^  =  0. 

Another  classic  invariant  is  the  energy  E  (per  kilogram)  defined  by 

£■  =  •  V  —  k/r 

An  easy  computation  shows  that  ^  =  0  and 

2E  =  ||cl|2-i?2  =  _(i  _  e2)A:2//i2  (2.1) 

2.2  Characterization  of  Orbits  with  Fixed  Energy  and  Eccen¬ 
tricity 

Our  nominal  strategy  is  to  work  with  orbits  of  fixed  energy  and  eccentricity.  Here  we  describe  these 
orbits  in  terms  of  velocity-circle  coordinates  and  examine  impulsive  orbit  transfers  that  preserve 
these  two  parameters. 
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2.2.  Characterization  of  Orbits  with  Fixed  Energy  and  Eccentricity 


If  we  fix  both  E  and  e  of  the  orbits,  from  the  definition  of  €  =  |  |c|  |/i?,  we  find 

2E  =  R^(€^  -  1) 

Thus, 

R  =  sj2El{e^  -  1)  and  ||c||  =  ci? 

show  that  both  ||c||  and  R  are  uniquely  determined.  Also  note  that  the  magnitude  of  the  angular 
momentum /i  =  A://?isdetennined. 

These  equations  clearly  identify  geometrically  the  set  of  orbits  of  fixed  energy  and  eccentricity. 
In  the  velocity  space,  the  four-parameter  set  of  orbits  is  obtained  as  follows: 

1.  Start  with  the  velocity  circle  of  any  one  such  orbit. 

2.  Apply  an  orthogonal  transformation  (three  parameters). 

3.  Rotate  the  circle  about  its  center  (one  parameter). 

The  orthogonal  transformations  alter  the  orbital-plane  parameters  (classically,  these  are  the  incli¬ 
nation,  longitude  of  the  ascending  node,  and  argument  of  perigee)  and  the  direction  of  the  satellite 
along  the  orbit.  The  rotation  about  the  center  of  the  velocity  circle  is  a  shift  of  the  true  anomaly  B. 

We  can  now  characterize  those  impulsive  transformations  that  preserve  both  energy  and  ec¬ 
centricity.  First,  because  the  transfer  is  impulsive,  the  position  vector  x  remains  fixed  (only  the 
velocity  v  changes).  But  now: 

1.  The  energy  is  fixed,  so  |  |v|  |  is  constant. 

2.  Because  h  is  constant,  ||x  x  v||  is  constant. 

The  set  of  v  restrictions  satisfying  1  and  2  lies  in  two  circles  in  planes  orthogonal  to  x.  The 
only  transformations  of  the  velocity  vector  satisfying  1  and  2  simultaneously  are  composed  of  the 
following  two  transformations: 

1 .  Rotate  the  velocity  vector  by  an  arbitrary  angle  5  about  the  x  vector  (the  direction  in  the 
plane  of  the  velocity  circle  through  0  and  orthogonal  to  v  —  c). 

2.  Reflect  v  with  respect  to  the  plane  perpendicular  to  x. 

Transformations  of  the  first  type  are  reasonable  and  easy  to  compute  mathematically.  They  can  be 
performed  practically  for  small  values  of  S.  Transformations  of  the  second  type  may  be  feasible  if 
the  angle  between  v  and  x  is  close  enough  to  90  degrees. 

In  summary,  we  can  explicitly  parameterize  the  velocity  circles  corresponding  to  orbits  of 
fixed  energy  and  eccentricity  as  the  product  space  of  the  three-dimensional  orthogonal  group  and 
the  circle.  Furthermore,  it  is  easy  to  compute  (using  rotations  and  reflections)  all  single-impulse 
AV  commands  that  preserve  the  same  two  parameters. 
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2.3.  Separation  Function 


2.2.1  More  General  Orbit  Transfers 

To  parameterize  the  more  general  set  of  single-impulse  orbit  transfers,  in  which  it  is  desired  to 
correct  the  total  energy  and  eccentricity  as  well  as  adjust  the  range  of  a  separation  function,  start 
by  finding  a  single  Vnew  having  the  desired  E  and  e  parameters  for  the  current  x.  To  do  this, 
we  need  to  solve  two  quadratic  equations  in  three  unknowns.  Geometrically,  Vnew  lies  in  the 
intersection  of  a  cone  and  a  sphere,  easily  found  by  analytic  methods: 

||v„ew|P  =  2(£;  + A:/||x||) 


I  |x  X  V new  1 1  —  ^netu  — 

To  obtain  a  more  general  set  of  transformations,take  any  solution  Vnew  to  these  equations  and  apply 
to  it  the  general  set  of  transformations  of  Vnew.  preserving  E  and  e  as  discussed  in  the  previous 
subsection,  to  obtain  the  more  general  set  of  transformations. 

2.3  Separation  Function 

A  separation  function  is  a  map  from  pairs  of  satellites  to  R^,  the  real  plane.  The  first  coordinate 
is  an  estimate  of  how  close  the  two  satellites  will  come  during  the  next  orbital  period,  the  second 
coordinate  is  an  estimate  of  how  far  apart  they  will  drift. 

A  good  separation  function  has  three  properties: 

1 .  The  first  value  vanishes  whenever  the  two  satellites  are  on  a  collision  course. 

2.  The  second  value  is  approximately  equal  to  the  maximum  separation  distance  during  the 
next  orbital  period. 

3.  It  is  analytically  tractable. 

The  significance  of  the  first  two  properties  should  be  clear.  Indeed,  the  most  obvious  selection 
for  the  separation  function  (called  the  obvious  function)  is  constmcted  as  follows:  simulate  the 
Kepler  orbits  of  the  two  satellites  over  the  next  orbital  period  and  compute  the  maximum  and  min¬ 
imum  distance  separations  during  that  period.  Unfortunately,  the  obvious  function  is  less  tractable 
than  we  want,  thus  violating  the  third  property.  From  a  theoretical  perspective,  the  difficulty  with 
the  obvious  function  is  the  transcendental  nature  of  the  true  anomaly  as  a  function  of  time  (except 
when  c  =  0,  discussed  below). 

As  a  partial  alternative  to  the  obvious  function,  we  are  experimenting  with  a  candidate  first 
coordinate  for  a  separation  function  called  the  cross-path  delay.  The  cross-path  delay  is  defined  as 
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2.3.  Separation  Function 


follows:  given  the  two  orbits,  compute  the  line  on  which  the  two  orbital  planes  intersect.  There  are 
two  pairs  of  points,  {xia,X2a)  and  (xib,X2b)  (subscript  1  corresponds  to  the  first  orbit,  subscript 
2  to  the  second),  near  to  each  other  on  this  line  of  intersection.  Note  that  if  the  centers  of  gravity 
of  the  two  satellites  are  going  to  collide  anywhere,  it  would  have  to  be  at  one  of  these  two  pairs  of 
points.  The  cross-path  delay  of  the  first  pair  of  points  is  the  difference  in  the  times  when  satellite  1 
will  be  at  xia  and  when  satellite  2  will  be  at  X2a-  The  cross-path  delay  of  the  second  pair  is  defined 
similarly,  and  the  overall  cross-path  delay  is  the  minimum  of  these  two  values. 

There  are  two  attractive  features  of  the  cross-path  delay  function.  First,  it  should  be  a  good 
measure  of  the  apparent  separation  (as  seen  from  the  Earth)  between  the  two  satellites  when  they 
approach  the  crossing  nodes.  Second,  the  cross-path  delay  is  a  readily  computable  function.  The 
approach  to  the  computation  is  as  follows:  note  the  current  true  anomalies,  OiandOi,  of  the  two 
satellites  and  compute  the  true  anomalies  corresponding  to  the  crossing  points  (two  on  each  orbit). 
It  is  possible  to  derive  an  analytic  expression  for  time  as  a  function  of  the  trae  anomalies.  Given 

d9  _  h{l  +  €Cos{9))^ 
di  ~  A2 

we  can  rewrite  this  as 

f.  A2  /•  d9 

~  h  I  {I  +  ecos(0))2 

and  the  right  side  of  this  equation  (2.2)  can  be  evaluated: 

r  d9  6sin(0)  2  \/l  -  e^tand) 

y  (1  -I-  ecos(0))2  (e2  -  i)(i-f  ecos(0))  (i  —  e^)!  1  -f  c 

The  cross-path  delay  can  be  computed  accurately  from  this  expression.  For  orbits  having  ec¬ 
centricities  around  0.0001  or  smaller,  there  are  approximate  solutions  to  the  above  integral  that  are 
far  easier  to  evaluate.  The  approximate  solution  to  first  order  in  e  is 

j  dt  =  ^(0  —  2esin(0))  ■ 
while  the  solution  to  second  order  in  c  is 

j  dt  =  '^{9-2esm{9)  -I-  3c2(^ -|- ?H^®)) 
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Section  3  Satellite  Cluster  Reconfiguration 


In  the  course  of  a  mission  such  as  that  outlined  in  the  TechSat21  program,  a  satellite  cluster  may 
have  to  undergo  repeated  reconfiguration. 

Repeated  reconfiguration  of  satellites  presents  the  problem  of  equalization  of  fuel  use  across 
the  cluster.  Naive  reconfiguration  algorithms  (e.g.,  each  satellite  has  a  fixed  station  that  it  takes 
at  reconfiguration  time)  may  produce  very  unequal  fuel  use,  resulting  in  a  suboptimal  total  clus¬ 
ter  lifetime.  We  wish  to  equalize  fuel  use  across  a  sequence  of  many  maneuvers  and  do  it  in  a 
distributed  fashion  requiring  minimal  computational  resources.  We  assume  that  for  n  satellites 
in  a  cluster,  there  are  n  known  orbital  stations  to  be  occupied  in  a  reconfiguration,  and  that  each 
satellite  can  calculate  its  cost  to  move  to  each  of  the  stations.  Our  problem  thus  becomes  assigning 
satellites  to  stations  so  as  to  maximize  cluster  lifetime. 


3.1  Market-Oriented  Programming 

Wellman  [6]  has  developed  an  auctioneering  model  based  on  the  economic  theory  of  general  equi¬ 
librium.  His  model  has  been  applied  to  problems  of  transportation,  computing  resource  allocation, 
and  product  design  (see  for  example  [5]). 

We  investigated  the  use  of  tnarket-oriented programming  techniques  to  select  the  stations  that 
satellites  occupy  during  reconfigurations.  In  our  study  problem,  we  have  a  set  of  n  homogeneous 
agents,  the  satellites,  each  of  which  must  “produce”  one  of  n  goods,  i.e.,  stations.  There  is  just 
one  resource,  fuel,  and  each  satellite  has  a  fixed,  remaining  allocation  of  that  resource;  they  cannot 
trade  fuel  among  themselves.  For  this  reason,  we  decided  to  begin  with  a  highly  simplified  auction 
algorithm.  In  this  algorithm,  each  satellite  simply  submits  a  fixed  bid  for  each  station.  This  bid 
is  based  on  a  cost  function  of  the  following  form  Bij  =  where  Bij  is  the  bid  submitted 

by  satellite  i  for  station  j,  Auj j  is  the  Au  required  for  satellite  i  to  reach  station  j,  and  r/j  is  the 
remaining  fuel  for  satellite  i. 

All  satellites  bids  are  then  submitted  to  an  auction  (which  might  be  mn  on  one  of  the  satellites  or 
at  a  ground  station).  This  auction  then  assigns  satellites  to  stations  using  the  following  algorithm: 

1 .  Determine  the  lowest  bidder  for  all  stations  yet  to  be  assigned. 
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3.2.  Choosing  a  cost  function 


2.  Determine  the  maximum  among  stations  of  the  lowest  bids. 

3.  Assign  the  station  corresponding  to  this  maximum  to  its  lowest  bidder. 

4.  Remove  the  assigned  station  and  satellite  from  the  auction. 

5.  If  not  all  stations  have  been  assigned,  go  to  step  1. 

3.2  Choosing  a  cost  function 

As  discussed  above,  the  cost  function  used  to  determine  bids  must  increase  as  the  Av  required 
increases  and  must  increase  as  the  fuel  remaining  decreases.  Therefore,  the  simplest  possible 
function  is  Bij  =  i.e.,  the  bid  is  simply  the  cost  as  a  fraction  of  fuel  remaining.  We  could 

also  introduce  an  exponential  dependence  on  fuel  remaining,  for  example,  Bij  =  or  an 

exponential  dependence  on  Av:  Bij  =  of  course,  there  are  numerous  other  variations 

possible.  We  tested  the  above  three  functions  to  determine  whether  using  an  exponential  depen¬ 
dence  on  one  of  the  two  main  parameters  increased  the  effectiveness  of  fuel  equalization. 


3.3  Test  problem 

As  an  initial  test  of  our  auction  algorithm,  we  used  the  following  reconfiguration  maneuver:  move 
from  a  parking  orbit  (i.e.,  all  satellites  in  the  same  circular  orbit.  Figure  3.1)  to  the  mission  con¬ 
figuration  (Figure  3.2).  The  parking  orbit  consisted  of  17  satellites  spaced  out  in  a  circular,  polar 
orbit  with  an  altitude  of  637  km  and  an  intersatellite  separation  of  44  m.  (To  gauge  the  effects  of 
different  perturbations  we  added  ,one  satellite  at  the  constellation’s  center). 

For  a  single  satellite,  then,  a  transfer  from  the  parking  orbit  to  one  of  the  elliptical  orbits  in  the 
“circle”  required  an  in-plane  bum  to  change  orbital  eccentricity  and  an  out-of-plane  bum  to  shift 
the  orbital  plane  very  slightly.  The  estimated  costs  for  these  maneuvers  ranged  from  4  cm/sec  to 
83  cm/sec.  Each  satellite  was  initially  assumed  to  have  sufficient  fuel  for  a  Av  of  100  m/s;  this  is 
consistent  with  a  fuel  mass  fraction  of  10%  and  a  specific  impulse  of  100  s. 

The  results  show  that  the  simple  auction  algorithm,  using  the  simplest  cost  function  Avjj  /rfj, 
produced  a  near-optimal  assignment  of  configurations  among  the  satellites.  The  satellites  were  able 
to  perform  242  reconfigurations  before  one  of  them  ran  out  of  fuel,  and  after  those  reconfigurations 
only  0.6%  of  the  total  fuel  was  left  unused.  By  comparison,  a  naive  algorithm  that  sets  bids  equal 
to  Av  regardless  of  fuel  remaining,  and  thus  always  produces  the  same  assignment  of  satellites  to 
stations,  resulted  in  a  satellite  running  out  of  fuel  after  141  reconfigurations,  with  42.3%  of  the 
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3.3.  Test  problem 


Figure  3.1:  Formation  before  reconfiguration  (distances  in  km) 

total  fuel  unused.  The  fuel  usage  histories  for  the  satellites  in  these  two  tests  are  plotted  in  Figures 
3.3  and  3.4. 

The  cost  function  exp(A'Uij) /r/j  still  ran  out  of  fuel  after  242  reconfigurations,  but  achieved 
a  slightly  more  uniform  fuel  distribution,  with  only  0.4%  of  the  fuel  left  unused.  The  function 
Aujj/  exp(r/i),  however,  did  very  poorly,  mnning  out  of  fuel  after  142  reconfigurations  with 
41.8%  of  the  fuel  left  unused.  A  slight  variation  on  this  function,  setting  the  bid  to  Avij  x 
exp(l  /rfi),  achieved  much  better  results.  It  ran  out  of  fuel  after  239  configurations  with  only 
1 .9%  of  the  fuel  still  unused.  The  fuel  usage  histories  for  these  tests  are  plotted  in  Figures  3.5, 3.6, 
and  3.7.  All  of  these  results  are  summarized  in  Table  3.1. 

Table  3.1 :  Summary  of  different  bidding  strategies 


Cost 

Function 

#of 

Reconfig. 

%  of  Fuel 
rentaining 

Avij/rfi 

242 

0.6 

Av 

141 

42.3 

exp(Ar;ij)/r/j 

242 

0.4 

Avij/  exp(r/i) 

142 

41.8 

Avij  X  exp(l/r/j) 

239 

1.9 
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Figure  3.2:  Formation  after  reconfiguration  (distances  in  km) 

3.4  Distributed  Auctioneering  for  Fault  Tolerance 

The  1999  TechSat21  program  developed  an  auctioneering  algorithm  to  optimize  the  station  keeping 
and  constellation  reconfiguration  decisions  and  to  maximize  the  lifetime  of  the  constellation.  The 
2000  auctioneering  work  focused  on  the  distribution  of  that  algorithm  to  provide  fault  tolerance 
and  load  balancing,  while  minimizing  communication  and  computation  overhead. 

This  study  covered  several  distinct  distribution  mechanisms.  The  best  mechanism  depends  on 
the  number  of  satellites,  the  level  and  type  of  fault  tolerance  desired,  as  well  as  the  available  com¬ 
munication  and  computing  bandwidth.  For  instance,  a  very  small  constellation  of  2  to  3  satellites 
would  probably  use  a  simple  distribution  scheme  with  limited  fault  coverage.  A  medium  sized 
constellation  of  4  to  36  satellites  might  use  a  scheme  that  is  Byzantine  resilient.  Extremely  large 
constellations  would  probably  apply  a  hierarchy  of  approaches,  due  to  the  n3  overhead  of  the 
Byzantine  approach. 

The  key  assumptions  in  this  effort  are  as  follows: 

1.  Relative  satellite  positions  are  available  (GPS) 

2.  Sufficiently  precise,  synchronized  time  is  available  (GPS). 

3.  Significant  on-board  processing  is  available. 

4.  Satellites  have  homogeneous  capabilities. 

5.  A  reliable  communications  infrastructure  is  in  place  (e.g.,  a  garbled  message  is  detectable). 
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FubI  um  hnlory.  oosi  =  (Ma-v  t  fuel  retnaininQ 


0  50  10O  150  200  -  250 

Nu(nt>*r<rf  r*oonAguralbm  p«r1onn«d 


Figure  3.3:  Fuel  usage  history  with  cost  function  Avijfrfi 

A  distribution  algorithm  is  considered  “safe”  if,  for  a  given  type  of  failure,  two  satellites  cannot 
attempt  to  assume  the  same  station  in  the  constellation.  The  failure  modes  that  were  considered 
include,  but  are  not  limited  to  the  following:  passive  satellite  failure,  communication  problems 
(detectable  by  the  communication  subsystem),  loss  of  GPS  (time  or  position  errors),  sensor  failures 
(i.e.,  errors  in  position  or  fuel  measurement),  byzantine  (asymmetric)  failures.  For  example,  if  an 
errant  satellite  A  sent  different  (but  self-consistent)  status  information  to  satellites  B  and  C,  then 
B  and  C  could  come  to  conflicting  solutions.  It  has  been  proven  that  systems  with  less  than  four 
contributors  and  two  exchanges  cannot  be  made  resilient  to  this  class  of  fault  (see  [2]). 

The  following  are  the  reconfiguration  approaches  analyzed  under  this  program: 

Master/Slave:  In  this  approach,  developed  under  the  1999  TechSat21  program,  a  master  col¬ 
lects  bids  from  everyone  and  distributes  the  results.  This  is  not  resilient  to  arbitrary  failures,  passive 
or  otherwise. 

Distributed  Minimized  Computation:  This  approach  distributes  the  proposed  master/slave  al¬ 
gorithm,  which  handles  arbitrary  fail  passive  behavior.  It  attempts  to  minimize  overall  computation 
cycles  by  having  each  satellite  receive  bids  from  all  other  satellites,  do  the  bid  selection  itself,  and 
then  act  on  the  result. 

Distributed  Minimized  Communication:  This  approach  distributes  the  proposed  master/slave 
algorithm,  which  handles  arbitrary  fail  passive  behavior.  It  attempts  to  minimize  communication 
overhead  by  having  each  satellite  broadcast  its  state.  Each  satellite  then  receives  these  states. 
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Figure  3.4:  Fuel  usage  history  with  cost  function  Av 

computes  bids  for  itself  and  the  other  satellites,  does  the  bid  selection,  and  then  acts  on  the  result. 

Distributed  Dependable:  This  approach  distributes  the  resulting  assignment  vector  from  a  dis¬ 
tributed  algorithm,  includes  total  cost  with  the  assignment  vector,  and  selects  the  vector  that  is  re¬ 
peated  most  often.  It  then  breaks  ties  by  selecting  the  lowest  cost  solution  and  then  breaks  ties  with 
that  by  selecting  the  one  provided  by  the  lowest  identification  number.  This  one  is  Byzantine-fault 
secure  for  formations  of  four  and  greater.  Byzantine  resilience  needs  k+1  rounds  of  information 
exchange  to  tolerate  k  faults.  This  approach  has  two  rounds  and  is  therefore  tolerant  to  just  one 
Byzantine  fault. 

Distributed  Byzantine  Resilient:  This  approach  (which  can  also  optimize  for  plume  avoidance) 
uses  the  distributed  dependable  algorithm,  but  distributes  the  result  after  each  assignment  and 
collects  new  bids  following  each  assignment.  This  enables  enhancements  for  handling  additional 
Byzantine  faults  and  also  allows  path  calculations  for  collision  and  plume  avoidance  to  be  factored 
into  cost.  At  each  step,  this  approach  only  involves  the  satellites  that  do  not  yet  have  station 
assignment,  which  provides  overall  lower  communication  costs. 

Another  approach  is  to  combine  instances  of  the  above,  for  example,  distributed  or  master/slave 
at  one  level  and  Byzantine  resilient  at  another.  This  is  potentially  attractive  since  the  Byzantine- 
resilient  algorithms  have  significant  nonlinear  overheads  so;  laige  clusters  might  be  unachievable 
without  a  hybrid.  For  instance,  the  Byzantine  solution  for  a  256-element  formation  is  126  times 
slower  than  the  proposed  minimally  communication-intensive  algorithm,  and  it  requires  1840  times 
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Figure  3.5:  Fuel  usage  history  with  cost  function  exp{Avij)/rfi 


the  communication  bandwidth. 

However,  these  hierarchical  approaches  suffer  from  additional  design  complexity  and  increased 
cost  for  smaller  constellations  over  the  more  straightforward  approach.  Issues  include  domain 
selection  and  arbitration.  Examples  of  such  hierarchical  approaches  include; 

1 .  Holding  an  election  for  a  master. 

2.  Holding  an  election  for  several  masters. 

3.  Voting  results  with  Byzantine  exchange. 

4.  Applying  a  Byzantine  approach  to  a  cluster  of  satellites,  then  applying  a  Byzantine  approach 
to  the  clusters  (this  has  unresolved  issues  with  cluster  selection). 

5.  Bidding  for  a  region  (stations  <C  satellites),  with  multiple  winners  per  region  making  up  a 
cluster,  then  bidding  for  spots  within  that  cluster. 

6.  Clustering  based  on  location  in  the  current  configuration  (some  locality  benefits,  but  has 
potential  boundary  effects.) 
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Section  4  Formation  Control 


4.1  Perceptive  Control  Theory 

Our  formation  control  layer  algorithms  are  based  on  perceptive  control  theory  [1].  This  theory 
has  been  successfully  applied  to  formation  control  of  multiple  autonomous  agents.  In  perceptive 
control  theory,  the  key  component  of  a  system,  such  as  a  formation  of  multiple  satellites,  is  defined 
to  be  the  perceptive  action  reference.  This  action  reference  is  computed  on-line  based  on  sensory 
measurements.  The  on-line  planner  of  the  system  generates  the  desired  state  values  for  the  system 
according  to  the  computed  action  reference.  In  addition,  the  action  reference  is  calculated  near 
or  at  the  same  rate  as  the  feedback  control.  In  other  words,  the  action  plan  is  adjusted  at  a  high 
rate,  which  enables  the  planner  to  handle  unexpected  or  uncertain  discrete/continuous  events  or 
planned  system  reconfigurations.  Furthermore,  unlike  traditional  methods  that  require  controller 
replanning  to  complete  the  task,  once  the  event  or  reconfiguration  is  over,  the  perceptive  control  in 
this  case  does  not  require  any  replanning. 

A  formation  of  multiple  spacecrafts  has  a  special  action  reference.  The  tasks  are  synchronized 
according  to  the  given  action  reference.  A  perceptive  frame  will  be  a  mathematical  abstraction  of 
these  action  references.  It  is  a  projection  from  the  position  p  eV  and  the  force  f  e  iF  to  a.  reference 
s  e  S,U  :  V  X  T  S.  Based  on  this  projection,  an  action  plan  of  a  multi-vehicle  formation  in 
the  perceptive  frame  is  parameterized  by  the  corresponding  action  reference  parameter.  Since  the 
action  reference  is  a  function  of  the  real-time  sensory  information,  the  desired  quantities  generated 
by  the  action  plan  are  also  directly  related  to  the  measured  data.  This  creates  a  mechanism  to 
provide  information  in  the  form  of  numerical  values  for  the  base  plan  using  the  measurements, 
and  the  information  is  interpreted  through  a  perceptive  frame  to  determine  the  control  input  value. 
Thus,  the  planning  becomes  a  closed-loop,  real-time  process.  The  plan  (desired  action)  can  be 
perceived  as  an  “abstract”  entity,  like  a  function.  It  only  has  a  “real”  value  when  measurements  are 
made  within  the  given  perceptive  frame. 
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4.2.  Application  of  Perceptive  Control  Theory  to  Cluster  Formation  Flying 


4.2  Application  of  Perceptive  Control  Theory  to  Cluster  For¬ 
mation  Flying 

Perceptive  control  theory  is  a  general  method  for  controller  design  of  formations  of  multiple  au¬ 
tonomous  agents.  A  simplified  schematic  of  the  control  architecture  is  shown  in  Figure  4.1. 


Formation  center 


Synthetic  time 


Figure  4.1:  Formation  control  layer  architecture 

To  develop  a  formation  control  scheme  for  multiple  spacecraft  vehicles,  the  perceptive  control 
method  will  be  generalized  to  a  broader  sensing-based  perceptive  context.  The  first  step  is  to 
find  an  appropriate  perceptive  action  reference  that  can  efficiently  represent  and  transfer  the  output 
measurements  of  each  vehicle  to  all  action  planners.  The  second  step  is  to  develop  a  computational 
scheme  for  determining  the  perceptive  reference  based  on  real-time  sensor  measurements.  The 
third  step  is  the  design  of  a  general-purpose  controller  with  respect  to  the  perceptive  reference 
frame. 

In  our  approach,  we  will  compute  a  perceptive  reference  variable — synthetic  time — ^that  is 
communicated  to  the  individual  vehicles,  which  will  use  this  variable  to  generate  local  guidance 
commands.  Given  the  current  state  of  the  set  of  vehicles  si  {t)  ,S2(t),...,  s„(t)  and  their  preplanned 
trajectory  indexed  in  a  parameter  r,  tTi(r),  <72(t),  . . . ,  cr„(r),  the  reference  variable  &  is  computed 
as  0  =  argmin.r(ll(5i(t), 52(f),..., Sn(f))  -  (<^i(’’),<72(r), . . . ,<r„(r))i|p).  It  can  be  shown  that 
for  norms  F  that  meet  a  certain  set  of  conditions,  the  system  represented  in  Figure  4.1  is  stable 
in  the  sense  of  Lyapunov  if  the  individual  components  are  stable.  Once  0  is  determined,  each 
individual  satellite  regulates  to  the  command  cri(ff).  The  nominal  trajectories  can  be  stored  in  each 
vehicle  relative  to  a  cluster  center.  In  this  way,  the  trajectory  is  separated  into  two  components, 
one  describing  the  overall  behavior  of  the  cluster  and  one  describing  its  shape. 


Formation  center 


Synthetic  time 


Formation  center 


Synthetic  rime 


Formation  center  I 


Synthetic  time 
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Reference  Projection 
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4.3.  Virtual  Leader  Reference  Projection 


4.3  Virtual  Leader  Reference  Projection 

For  our  example,  we  use  the  position  of  the  “center”  of  the  constellation  as  the  reference  variable. 
Sincp  Jie  trajectories  we  are  considering  are  symmetric,  we  computed  the  center  by  taking  the 
average  of  the  current  positions.  Let  Pi  be  the  absolute  position  of  the  satellite  at  time  t,  and 
lei  Pni{B)  be  a  parameterization  of  the  nominal  trajectory  of  the  satellite.  Then,  the  measured 
center  at  time  t  is  given  by  (7(f)  =  ^  D"  Pu  and  a  parameterization  of  the  nominal  trajectory 
for  the  constellation  center  is  CniO)  =  ^  The  reference  projection  is  then  obtained 

by  finding  the  point  in  the  nominal  trajectory  for  the  constellation  center  closest  to  the  current 
center.  We  used  as  distance  the  geometric  distance  (i.e.,  only  the  position  states  were  used  in  the 
projection).  Figure  4.2  shows  how  this  is  done  for  the  case  of  a  constellation  of  two  satellites. 
The  reference  variable  will  be  the  time  index  9  that  minimizes  the  distance  between  the  nominal 
center  at  time  0  and  the  current  position  of  the  center.  We  will  denote  this  projection  by  11.  9  = 
n(C')  =  argmin  (||(7(f)  —  (7„(0)||) .  This  approach  has  the  advantage  of  not  requiring  the  setup 
of  a  master/slave  architecture.  We  are  thus  able  to  have  a  completely  symmetric  formation  control 
layer.  The  benefits  of  this  are  simpler  architecture  and  more  fault  tolerance. 


Nominal  Poeition  at  synthetic  time 

Figure  4.2:  Reference  projection  using  constellation  center 


4.4  Use  of  Relative  Position  Sensing 

Relative  position  between  satellites  can  be  obtained  with  much  higher  accuracy  than  absolute  po¬ 
sition.  Since  we  are  tracking  a  nominal  trajectory  given  in  absolute  position  for  each  satellite,  it 
may  seem  like  we  will  suffer  high  errors.  However,  by  expressing  the  nominal  trajectories  for  each 
satellite  relative  to  the  constellation  center,  and  by  computing  the  constellation  center  using  both 
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4.4.  Use  of  Relative  Position  Sensing 


relative  and  absolute  position,  we  can  easily  work  around  this  problem.  In  this  case,  each  satellite 
will  compute  a  reference  projection.  As  we  shall  see,  this  will  not  significantly  affect  the  shape  of 
the  formation.  If  Pi  is  the  absolute  position  of  the  satellite,  the  relative  position  between  the 
and  satellites  will  be  given  L-y  Pij  =  Pj  —  Pi.  Let  Cj  be  the  error  in  absolute  position  for  the 
satellite.  Then  the  satellite  can  compute  the  measured  center  of  the  constellation  as; 

Cmi  =  -(n{Pi  +  ei)  +  j2Pij 

”  \  i=2 

The  satellite  now  computes  its  position  command  P?  by  summing  the  desired  relative  positions 
from  nominal  center  to  the  measured  position.  Knowing  past  inputs  and  outputs  allows  us  to 
estimate  current  state  and  model  parameters  position  of  the  center  Tf  =  n  (II  (C  +  ej))+C+ej.  We 
can  approximate  this  expression  by  linearizing  the  projection,  P^  «  r;  (n(C))  +  Vri(n(C))ei  + 
C  +  Cj.  The  measured  tracking  error  e  will  now  be  e  =  P/  -  Pi^  «  e  +  C{t)  —  C  (n(C))  + 
Vfj  (n  (C))  Ci,  where  e  is  the  actual  tracking  error.  Forcing  the  measured  tracking  error  to  be 
zero,  e  =  0  =»  e  =  —  {C{t)  —  C  (n(C))  +  Vrj  (E  (C))  Ci) .  The  real  tracking  error  thus  has 
two  components:  one  that  doesn’t  depend  on  the  satellite  and  thus  doesn’t  affect  the  constellation 
shape,  and  one  that  depends  on  the  measurement  error  and  the  gradient  of  the  nominal  trajectory. 
Since  this  will  be  a  small  number,  it  will  dampen  the  effect  of  the  measurement  errors. 
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Section  5  Tradeoff  Studies 


In  this  section  we  explore  the  effect  of  different  factors  on  the  overall  fuel  usage  of  the  constellation 
while  maintaining  formation.  To  prepare  for  the  demonstration  flight  of  the  Techsat21  mission  in 
2003,  it  is  necessary  to  establish  which  are  the  most  fuel-efficient  strategies  to  operate  the  cluster 
while  in  formation.  It  is  also  important  to  determine  under  which  strategy  the  demonstrator  mission 
will  last  up  to  a  year.  Of  that  length  of  time,  only  a  small  fraction  will  be  used  to  test  formation 
flying  capabilities.  The  total  fuel  budget  for  each  satellite  in  the  demonstrator  mission  is  equivalent 
to  a  total  AV  of  65m/s. 

A  key  metric  for  the  control  strategy  is  thus  the  constellation  endurance,  measured  by  the 
length  of  time  that  the  cluster  can  be  maintained  in  a  given  formation.  All  test  results  are  for  a 
reference  orbit  having  a  400  km  altitude,  30  deg  inclination,  and  a  2  mm/s  relative  velocity  error, 
with  four  satellites  per  cluster.  We  tested  two  formations  centered  around  this  reference  orbit;  an 
out-of-plane  formation,  where  the  projection  of  the  satellites  on  a  horizontal  plane  centered  at  the 
reference  orbit  lie  in  a  circle  of  1km  diameter  (Figure  5.1),  and  an  in-plane  formation,  where  all 
satellites  lie  in  an  ellipse  centered  on  the  reference  orbit  and  with  semi-major  axes  of  1  km  along 
track  and  2  km  in  the  radial  direction.  These  are  challenging  formations,  likely  to  be  more  stringent 
than  what  would  be  needed  for  mission  purposes. 


5.1  Effects  of  Correction  Strateg3-^ 

An  important  factor  on  the  total  fuel  usage  of  the  constellation  is  the  rate  at  which  corrections  to 
the  orbits  are  planned  and  executed.  We  analyzed  the  effects  of  two  parameters:  the  time  between 
orbit  corrections  and  the  time  allowed  for  the  correction.  In  the  following  figures,  we  present  the 
results  obtained  with  two  of  the  strategies  analyzed.  One  is  a  correction  every  two  orbits,  allowing 
for  a  whole  orbit  to  finish  the  correction.  The  other  is  for  a  correction  every  orbit,  allowing  for 
a  half  orbit  to  finish  the  correction.  Figures  5.4  and  5.7  show  the  fuel  usage  for  each  of  the  four 
satellites.  Figures  5.2,  5.3,  5.5,  and  5.6  show  the  relative  distance  between  two  of  the  satellites  in 
the  constellation.  To  meet  the  performance  requirements  of  the  Techsat21  mission,  this  distance 
has  to  deviate  less  than  10%  from  the  nominal.  The  nominal  distance  for  Figures  5.2  and  5.5  is 
1.4  km.  For  Figures  5.3  and  5.6  it  is  2  km.  As  a  reference  we  can  compare  these  values  with 
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5.2.  Effects  of  Restricting  Bums  to  a  Short  Interval 


Satellite  2 


Figure  5.1:  Nominal  location  of  satellites  and  relative  distances 


those  corresponding  to  an  in-plane  formation,  which  minimizes  the  effects  of  J2  and  higher  order 
perturbations.  These  results  are  shown  in  Figures  5.8  and  5.9. 


5.2  Effects  of  Restricting  Burns  to  a  Short  Interval 

The  total  fuel  usage  of  the  constellation  will  also  be  influenced  by  operational  constraints  that  only 
allow  usage  of  the  engines  at  specific  times.  For  thermal  reasons,  the  engines  can  only  be  used  for 
a  limited  duration  (5  to  10  noin)  and  then  require  up  to  a  five  hour  cooldown. 

This  results  in  a  loss  of  flexibility  for  the  linear  program,  which  leads  to  lower  system  perfor¬ 
mance.  To  gauge  the  impact  of  this  effect,  we  ran  two  test  cases.  For  these  tests,  we  used  the 
control  strategy  that  gave  the  best  performance  tradeoff  in  the  unrestricted-engine-use  case. 

We  ran  these  tests  for  the  two  formations  discussed  earlier.  The  results  for  the  more  challenging 
out-of-plane  formation  are  shown  in  Figures  5.10,  5.11,  and  5.12.  These  are  to  be  compared  to 
Figures  5.2,  5.3,  and  5.4,  respectively.  The  results  for  the  in-plane  formation  are  shown  in  Figures 
5.13  and  5.14. 
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OietancA  between  saletlites  1  and  2  (loti) 


5.2.  Effects  of  Restricting  Bums  to  a  Short  Interval 


Figure  5.2:  Distance  between  satellites  1  and  Figure  5.3:  Distance  between  satellites  1  and 
2:  unrestricted  bums,  correction  every  other  3:  unrestricted  bums,  correction  every  other 
orbit  orbit 


0.09r-- . — 1 - 1 - 1 - 1 - r 


Figure  5.4:  DeltaV:  unrestricted  bums,  correction  every  other  orbit 
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Figure  5.5:  Distance  between  satellites  1  and  Figure  5.6:  Distance  between  satellites  1  and 
2:  unrestricted  bums,  correction  every  orbit  3:  unrestricted  bums,  correction  every  orbit 


5.3.  Effects  of  Relative  Velocity  Sensing  Error 


5.3  Effects  of  Relative  Velocity  Sensing  Error 

Another  important  factor  in  fuel  consumption  is  the  error  in  relative  velocity  sensing.  To  quantify 
this  effect,  we  ran  simulations  in  which  the  planner  had  access  to  the  exact  orbital  model  in  order 
to  exclude  effects  of  mismatched  dynamics.  The  baseline  relative  velocity  error  was' taken  to  be 
6mm/s.  This  is  an  optimistic  (though  possible)  value  for  the  error,  given  the  technologies  proposed 
for  the  demonstrator  mission. 

Figures  5.15,  5.16,  and  5. 17  represent  the  fuel  usage  of  all  four  satellites  for  cases  with  baseline 
noise,  double  baseline  noise,  and  half  baseline  noise.  The  total  fuel  usage  grows  linearly  with 
the  measurement  noise.  For  baseline  noise,  note  also  that  the  effect  of  dynamics  mismatch  in 
the  trajectory  planning  algorithms  (as  shown  in  the  previous  section)  is  similar  in  magnitude  to 
the  effect  of  noise.  Although  the  former  can  be  reduced  by  using  more  sophisticated  nominal 
trajectories,  the  latter  is  a  fundamental  limitation.  The  results  are  summarized  in  Table  5.1. 

Table  5.1 :  Effect  of  velocity  sensing  error 


Relative  Velocity 

Fuel  Usage 

Endurance 

Error  (mm/s) 

per  Orbit  (mm/s) 

(days) 

3 

9.6 

420 

6 

15.6 

260 

12 

25.9 

157 

5.4  Summary  of  Constellation  Endurance  for  the  Demonstra¬ 
tor  Mission 

The  purpose  of  this  study  was  to  analyze  the  feasibility  of  the  demonstrator  mission  objectives 
given  the  operational  constraints.  Since,  in  the  demonstrator  mission  for  Techsat21  the  large  out- 
of-plane  formations  are  only  to  be  flown  for  a  few  days,  the  results  presented  in  Table  5.2  are  very 
encouraging.  However,  when  these  simulation  experiments  were  run  we  did  not  have  the  final  set 
of  operation  guidelines  for  the  constellation  (which  are  somewhat  still  evolving).  In  particular,  we 
did  not  consider  the  effect  of  doing  bums  along-track  and  across-track  sequentially.  The  effect  of 
having  to  use  the  engine  for  short  intervals  was  also  likely  underestimated.  We  did  not,  however, 
have  enough  information  from  the  engine  supplier  to  fully  understand  these  usage  restrictions. 


Use  or  disclosure  of  data  contained  on  this  sheet  is  subject  to 
the  restrictions  on  the  dtle  page  of  this  document 


28 


5.4.  Summary  of  Constellation  Endurance  for  the  Demonstrator  Mission 


Figure  5.15:  Fuel  usage  with  relative  velocity  error  of  6  mm/s 


Table  5.2:  Summary  of  performance  with  high-fidelity  gravity  model 


Formation 

Bums 

AV  per  Orbit 
(mm/s) 

Endurance 

(Days) 

In-Plane 

Unrestricted 

8.5 

408 

Two  10-min 

20.8 

167 

Out-of-Plane 

Unrestricted 

21.8 

159 

Two  10-min 

38.2 

91 
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OotiaV  {m/9) 


5.4.  Summary  of  Constellation  Endurance  for  the  Demonstrator  Mission 


Figure  5. 16;  Fuel  usage  with  relative  velocity  Figure  5.17:  Fuel  usage  with  relative  velocity 
error  of  1 2  mm/s  error  of  3  mm/s 


Section  6  Trajectory  Optimization 

In  this  section  we  analyze  how  trajectories  that  meet  the  mission  constraints  can  be  generated. 
We  follow  an  approach  very  similar  to  that  presented  in  [4].  However,  we  used  a  different  set  of 
constraints  that  better  match  the  objectives  of  the  Techsat21  mission  and  our  guidance  laws. 

6.1  Trajectory  Parameterization 

lo  optimize  the  trajectory,  it  was  divided  into  n  intervals.  Within  each  time  interval,  each  of  the 
coordinates  of  the  orbit  was  fit  with  a  third-degree  polynomial. 

To  compute  initial  values  for  the  coefficients,  we  began  with  Euler  orbits  for  all  the  satellites, 
satisfying  the  Hill  conditions  for  periodic  relative  motion.  Given  the  initial  and  final  position  and 
velocity  for  each  interval,  it  is  possible  to  determine  the  initial  coefficients  of  the  polynomial. 
Assume  the  initial  and  final  positions  of  the  satellite  are  represented  by 

R0-C0  +  Citi  4-  02^1  -f  Czt\ 

Rf  —  Co  +  Cit2  +  Citl  -f  C^tl 

Taking  the  derivative  of  tne  position  polynomials  yields  equations  for  the  initial  and  final  ve¬ 
locities  of  the  satellites,  which  can  be  represented  by 

Vg  =  0  Cl  2c2fi  +  3c3fj 
V'jf  =  0  -|-  Cl  -|-  202^2  "I"  3C3f2 

If  Ro,Rf,Vo,Vf,  and  the  initial  and  final  times  for  the  interval  (fi  and  <2)  are  known,  it  is 
possible  to  solve  for  the  four  coefficients  of  the  polynomial  by  constructing  a  matrix  such  that 

-  1  ti  tl  f?  1  [Co  1  [  ■ 

1  #2  tl  ^2  _  Rf 

0  1  2ti  3ff  ^  C2  “  K 

.  0  1  2f2  3fi  J  [  C3  J  L  . 
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6.2.  Constraints 


If  this  multiplication  is  repeated  for  each  time  interval  in  the  orbit,  the  initial  coefficients  for 
the  entire  orbit  can  be  calculated. 

For  the  optimization  routine  to  perform  correctly,  it  was  necessary  to  scale  the  parameters 
describing  the  trajectories  so  that  they  all  had  values  in  similar  ranges  and  so  that  they  all  had 
similar'  effects  on  the  cost  function. 

To  achieve  this  we  introduced  a  set  of  new  variables  denoted  do,di,d2,  and  d^.  We  chose  these 
variables  such  that 

dx  dx  dx  dx 

ddo  ~  ddi  dd2  ^  ddz 

One  way  of  achieving  the  correct  scaling  is  through  the  following  definitions  (repeated  for  each 
coordinate,  for  each  satellite,  and  for  each  time  interval): 

Co  =  CQi  +  do 
Cl  =  Cli  + 

C2  =  C2i  + 

C3  =  Csi  + 

where  <2  is  the  length  of  each  interval. 

6.2  Constraints 

We  had  two  sets  of  constraints.  One  set  imposed 
imposed  operational  requirements  on  the  resulting  orbits  for  all  satellites. 

The  first  set  of  constraints  required  that  the  position  and  velocity  of  each  satellite  at  the  end  of 
one  interval  be  within  a  certain  epsilon  of  the  position  and  velocity  of  the  satellite  at  the  beginning 
of  the  following  interval. 


t2 

^2 

*2 


continuity  of  the  orbits,  while  the  other  set 


\R1-R2f  <0an(i  |Vi-V2|^<0 

The  other  set  of  constraints  assured  that  the  constellation  maintained  a  functional  shape.  For 
this  experiment,  we  used  a  constellation  of  four  satellites  such  that  their  projections  on  a  horizontal 
plane  centered  at  the  reference  lie  on  a  fixed  circle  and  form  at  every  time  a  square.  For  a  square 
formation  to  be  “rigid,”  it  is  necessary  to  fix  five  of  the  six  distances  between  the  vertices.  Thus, 
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6.2.  Constraints 


the  projection  of  the  formation  will  remain  in  the  same  approximate  shape  by  imposing  that  five  of 
the  six  dimensions  remain  within  a  certain  tolerance  of  their  nominal  values.  For  our  experiments, 
we  chose  the  five  distances  shown  in  Figure  6.1,  and  we  allowed  a  10%  deviation  from  the  nominal 
distance. 


1 


3 

Figure  6.1:  Constraints  necessary  to  fix  the  cluster  shape 
For  a  1km  orbit,  the  constraint  on  each  of  the  diagonal  distances  was 

1.8  <  d  <  2.2 

and  the  constraint  corresponding  to  the  edges  was 

(v/2  -  .l>/2)  <  d  <  (\/2  +  .l\/2) 

To  compute  these  distances,  first  we  determined  the  center  of  the  constellation  by  averaging  the 
position  of  all  four  satellites.  Then  we  projected  all  four  satellites  into  the  horizontal  plane  (normal 
to  the  radial  direction)  through  the  constellation  center.  For  a  half-orbit  planning  horizon,  with  180 
second  intervals,  there  were  161  constraints  ( n-1  constraints  for  each  of  the  five  shape  constraints 
plus  n-2  position  and  velocity  constraints  for  each  satellite). 
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6.3  Objective  Function 

Our  objective  was  to  minimize  the  amount  of  fuel  needed  to  maintain  a  formation  that  meets  the 
mission  requirements.  To  do  this,  we  used  the  “inversion”  process  described  in  [4]  to  compute  the 
fuel  necessary  for  a  given  foimation.  We  measured  fuel  in  units  of  AV.  To  further  simplify  the 
process,  we  computed  the  total  AV  in  one  time  interval  as  the  acceleration  due  to  the  engines  in 
the  middle  of  the  interval  times.  The  accelerations  were  measured  for  a  length  of  one  time  interval. 
For  each  of  the  satellites  and  for  each  of  the  coordinates,  we  had 


dn 


t/  =  [  0  0  2  6t]x 
We  compute  the  acceleration  due  to  gravity  (a^)  with  the  following  formulas; 
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Finally,  we  computed  our  objective  function  as  the  sum  of  the  squares  of  the  A  V’s  correspond¬ 
ing  to  all  satellites,  all  directions,  and  all  time  intervals. 


6.4  Optimization 

The  optimization  routine  used  in  this  project  was  fmincon  from  the  Matlab  optimization  package. 
Fmincon  uses  a  sequential  quadratic  programming  (SQP)  algorithm.  The  number  of  variables 
depends  on  the  time  interval  being  used.  We  ran  test  cases  for  different  planning  horizons.  For  a 
planning  horizon  of  a  half  orbit  cut  into  180  second  intervals,  we  had  642  decision  variables  (  4 
satellites  times  4  coefficients  for  each  direction  times  3  directions  (x,y,z)  times  (n-1)  where  n  is  the 
number  of  time  intervals.) 

The  results  obtained  for  the  half-orbit  planning  horizon  are  shown  in  Figures  6.2,  6.3,  6.4,  and 
6.5.  Note  that  the  maximum  AV  observed  for  one  orbit  is  0.9  mm/s.  This  corresponds  to  4500 
days  of  continuous  operations  in  this  1km  formation.  Of  course,  this  is  the  nominal  AV.  The 
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actual  fuel  usage  will  be  higher  since  it  is  significantly  influenced  by  the  effects  discussed  in  the 
previous  section.  Note  also  that  to  maintain  the  Hill  periodical  orbits  in  the  presence  of  J2  requires 
a  AV  of  40  mm/s  per  satellite. 
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Figure  6.2:  Fuel  optimal  trajectory  shape 


Figure  6.3;  Fuel  optimal  trajectory  distance  between  satellites  1  and  2 

These  results  confirm  the  feasibility  of  the  inversion  approach  to  finding  optimal  trajectories 
presented  in  [4],  even  in  the  presence  of  the  fairly  large  number  of  constraints  required  to  impose 
all  the  operational  requirements.  The  computation  time  for  the  problem  mn  was  several  hours. 
(However,  we  computed  the  constraint  gradients  numerically,  which  significantly  slowed  down 
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Figure  6.4:  Fuel  optimal  trajectory  distance  between  satellites  1  and  3 

the  computation.)  Also,  as  for  every  optimization-based  algorithm,  guaranteeing  convergence  to  a 
valid  answer  systematically  is  difficult. 
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Figure  6.5:  Fuel  optimal  trajectory  ilV  for  all  satellites 
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